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ABSTRACT 

We show how to estimate the enclosed mass from the observed motions of an ensemble of 
test particles. Traditionally, this problem has been attacked through virial or projected mass 
estimators. Here, we examine and extend these systematically, and show how to construct an 
optimal estimator for any given assumption as to the potential. The estimators do not explicitly 
depend on any properties of the density of the test objects, which is desirable as in practice 
such information is dominated by selection effects. As particular examples, we also develop 
estimators tailored for the problem of estimating the mass of the Hernquist or NFW dark 
matter haloes from the projected positions and velocities of stars. 
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1 INTRODUCTION 

Here, we consider the general problem of estimating the enclosed 
mass (or equivalently the gravitational potential) from kinematical 
data on tracer populations. In other words, suppose there are N test 
particles moving in a gravitational potential <p{r) generated by a 
mass density p{r). The data available to us are the instantaneous 
positions r, and velocities v, of the test bodies. However, it is only 
rarely that the full phase space information is available and often 
only components of position and velocity along the line-of-sight 
are measured. From these data, we wish to estimate the underlying 
gravitational potential or mass by a robust and unbiased statistical 
method. 

This problem has many applications in modern astrophysics 
- including estimating the mass of the Milk y Way and M31 
from the kinematics of dis t ant satellites l Little & Tremaindl 19871 : 
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estimating the 



IWilkinson & Evans! 1 19991 : IWatkins et al" 
mass of the haloes of dwarf galaxies from the stellar v elocities 
dStrigari et al]|2008t IWalker et alj|2009t IWolf et al.|p201oh . and es- 
timati ng the mass of galaxy groups and cl usters from their mem- 
bers dHeisler et aflll985l : ITuIIv et alj|2006l) . In fact, the kinemati- 
cal properties of tracer populations are one of the richest sources 
of data on the distribution of dark matter in galaxies and clusters. 
Therefore, it is important to extract as much information from the 
data as we possibly can. 

Given the significance of the problem, there has been surpris- 
ingly little effort on d eveloping the systematic theory of mass esti- 
mators. Early work dLimber & Mathewsl fl 960) exploited the virial 
theorem to obtain 
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for the mass enclosed by N test particles with line-of-sight veloc- 



ities V;,i and projected positions Ri o f each particle. The problem 
with this method was pointed out by Bahcall & Tremaine] dl98ll) . 
namely that the virial mass estimator is both biased and inefficient. 
These latter authors introduced the alternative projected mass esti- 
mator 

G \ e I GN A- 1 

i 

where C is a constant determined by the host potential and the ec- 
centricity of the orbits. 

There has also b een substantial previous work done on scale- 
free mass estima tors jwhitdl 198ll: iKulessa & Lvnden-Belllll992l: 
lEvans et alj|2003l) . Recently. IWatkins. Evans & Anl d2010h formal- 
ized and expanded the ideas from these previous papers, presenting 
a variety of mass estimators tailored to scale-free potentials and 
densities. These work by taking weighted averages of the combi- 
nations of velocities and positions that remain invariant under sim- 
ilarity transformations. None the less, it is also important to devise 
mass estimators that are optimized for more realistic and specific 
astrophysical potentials. For example, there are cosmological argu- 
ments that the dark halo densi ty is cusped like r -1 at small rad ii and 
falls off like r~ 3 at large radii dNavarro. Frenk & Whitell 19951) . It is 
natural to look for mass estimators that build upon this assumption 
at the very start. 

In this paper, we show how to find such mass estimators tai- 
lored for any given potential. We also find that our new estimators 
do not depend on the number density of the tracers at all. This is 
a real advantage - e.g., in the case of the Milky Way, the variation 
of the number density of the known satellite galaxies and globular 
clusters is dominated by the selection effects and the true number 
density can only be guessed at. 

This paper is arranged as follows. In Sect.|2] we develop some 
general theories on mass estimators, showing how to construct one 
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suitable for given potential. In Sect. [3] we give a few specific exam- 
ples for the cases of astrophysical interest, whilst Sect. [4] sketches 
the extension to projected data. The self-consistent case, when the 
potential and the density of the tracers are related through the Pois- 
son equation, is dealt with in Sect. [5] Finally, in Sect. [6] we provide 
a discussion and conclusions. Some applications of our estima- 
tors to the widely-u sed cosmological halo mode l of iNavarro et al.l 
d 19951) are found in lEvans. An & Deasor] d2010h . which should be 
considered as a companion to the current paper. 



2 THE THEORY OF MASS ESTrMATORS 

2.1 Jeans' Equation and the Virial Theorem 

Suppose our tracer population has a number density v(r) and a ra- 
dial velocity dispersion cr 2 (r) and is moving in a spherical dark halo 
potential (p(r), which by Newton's Theorem satisfies 

d(p _ GM(r) 
dr r 2 



(1) 



where M(r) is the enclosed halo mass within radius r. These quanti- 
ties are related to one another through the spherical Jeans equation 
that reads 
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where 
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is the so-called Binney anisotropy parameter for the spherical sys- 
tem. The typical application of the Jeans equations involves deriv- 
ing the potential and the dark halo mass profile from the observed 
behaviour of the tracer density and velocity dispersions ('Jeans 
modelling'). If the observations are composed of discrete sample 
datapoints, this is subject to the uncertainties related to the binning 
and requires large number of datapoints to extract any meaningful 
information. An alternative when only moderate number of data- 
points are available is to consider the system as whole such as uti- 
lizing the virial theorem. The relation between these two is most 
obvious in a spherical system, for which integrating the spherical 
Jeans equation essentially results in the scalar virial theorem. 

In order to see this, we start by noting that the spherical Jeans 
equation reduces to an exact differential form 

Id , GM 

-^ T (Qv<r 2 r ) = -v— (2) 
Q dr r 1 

by means of the integrating factor Q = Q(r) satisfying 
dlng(r) = 2£r) 
dr r 

Next we find that the relation between the local three-dimensional 
velocity dispersion and the radial velocity dispersion is given by 

2 dln(r 3 e-') 
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The three-dimensional velocity dispersion of the tracers within the 
sphere of the radius of r out is thus given by 



<v 2 > = drr 2 vcr 2 
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where 
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is the total number of the tracers between an inner r m and outer r out 
radius. Integrating by parts, setting r in = and also using equation 
(0 results in 
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given that <x 2 is not divergent as r — > 0. Here the last integral actu- 
ally defines the total potential energy of the tracers within the same 
sphere, i.e., 

Im , 2 GM /GM\ \W\ 

\W\ = An drrv ; ( ) = - — -. 
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Hence, equation ((5) reduces to 



<v 2 > 
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where 



v(r O ut)0" 2 (rout) 
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Note that v out is the mean number density of the tracers in the sphere 
of the radius of r out . 

Equation l[6ll is in fact equivalent to the statement of the scalar 
virial theorem for a pressure-supported spherical system as \{v 2 ) 
is basically the kinetic energy per tracer particle associated with 
the random motion. The presence of the boundary term (that is, 
the surface term, 3<p 2 ) is due to the hard cut-off at r = r ollt , which 
can correspond to the situation when the tracers are confined within 
the spherical radius of r out through the external pressure. However, 
it is usual to drop the boundary term if the system as a whole is 
considered. 



2.2 Tracer Mass Estimators 

The virial theorem (eq. |6) is traditionally used to estimate the to- 
tal mass of the system. The integral mean value theorem indi- 
cates that there exists a kind of 'mean radius' r within the inter- 
val bounded by the outer cut-off r ouC (i.e., < f < r out ) such that 
{GM/r) = GM(j)jf. Therefore, if one ignores the boundary term, 
one can relate the mass within the 'mean radius' to the total veloc- 
ity dispersion of the tracers, M(r) = f<v 2 >/G, which may also be 
suspected from a simple dimensional analysis. If the distribution of 
the gravitating mass is known or assumed, the definition of f can be 
made precise and furthermore M(r) can be scaled to provide the es- 
timate of M(r out ). That is to say, we can relate the total mass of the 
system to the certain integrals of kinetic properties of the tracers. 

However, this approach suffers from drawbacks related to the 
fact that the mass estimate depe nds on two separate averages (see 
e.g.. iBahcall & Tremainelll98ll) . This difficulty is partially over- 
come by the use of 'mass estimator', that is, an average of partic- 
ular combinations of kinetic properties of the tracers that directly 
relates to the total mass rather than to the potential energy, as does 
the virial estimator. However, derivation of the proper form of the 
mass estimator requires some analysis of the dynamics of the sys- 
tem. Here we still consider the simplest case of the spherical system 
traced by a non-rotating relaxed population in equilibrium. 

First we note that the virial theorem (with the boundary term 
dropped) indicates that (v 2 )/(v 2 ) = 1 where v c = (GM/r) 1 ' 2 is the 
circular speed of the potential. From this, one may naively expect 
that <v 2 /v 2 ) ss 1, but the distributed mass and tracers only make this 
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approximately so. However, with a proper weighting f(r), we can 
actually show that there exists a relation (fv 2 /v 2 ) = 1, which will 
be subsequently used to derive a proper mass estimator. 

Let us assume for the moment that the spherical dark halo 
profile M(r) is known. Then we can show that the proper weighting 
function is given by 



f(r) = 4 - 2/3(r) • 



din M(r) 



dlnr 

_ dln(r 4 g-'M-') _ QM d / r 4 \ 
~ dlnT " ~Ar\QM)' 



(7) 



Then the 'weighted' average of the tracer radial velocities v 2 in a 
spherical system given by 



lfv 2 r \ An f' 01 " fv<r 2 r r 3 
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Integrating by part leads to 
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where we have also used equation (2). 

F or any physical system, we apply the result of I An & Evansl 
d2009h to find that lim r _ r 4 vcr 2 /M = 00 Hence, setting r m = 
leads to the vanishing inner boundary term, and thus 



IfK 

I «2 



1 + 3- 



(10) 



This relation can be rearranged to yield the total dark halo mass 
M BUt = M(r out ) within radius r mt , 



GM, 



3r outT , 



where 



fi(r) = (0 < r < r out ) 



(11) 



(12) 



is the normalized dark halo mass profile function. Note that 
dln/j/dlnr = dlnM/dlnr, and thus f(r) can be evaluated if fi(r) 
is known without any reference to M oul . Consequently, if we are 
willing to assume the functional form of the halo mass profile jl(r), 
which is normalized to be /2(r out ) = 1, in a spherical region of inter- 
est, r < r out , then the total halo mass M out within the same spherical 
region can be estimated through a particular average of kinematic 



1 Naively, we have lim r _>o r 3 v = for any physical y since otherwise 
there would be an infinite mass concentration of the tracers at the cen- 



tre whereas lim r _>o r<r 
careful examination of 



/M is typically fin ite 1 An & Evans! [2009). More 



An & Evanj 1 2009) indicates that even for the ex 



ceptional case that lim r _>o rtr 2 /M diverges, the boundary term still van- 
ishes a s lim r _,o r 3 v = is always dominant. That is to say, we infer 
from I An & Evansl (2009) that, given M(r) behaving as ~ r'~" as r — > 0, 
a > 2f3 - 3 is the sufficient condition for this. However, if v ~ r~ y as r — » , 
we have 2/3<y<Q' + 2- the first inequality is due to lAn & Evansl J2006l) 
and the second to the fact that the tracers cannot be cusped steeper than the 
dark halo - and so it is met. 



properties of the tracers, which is in practice inferred from the cor- 
responding discrete sample mean, i.e., 



1 y nj 
GN ^ n 



n/(r,) i 



(13) 



possibly further adjusted by the boundary term (3r out G~ l g 2 ) if nec- 
essary. 

Here, the presence of the boundary term is again related to 
the external pressure support of the tracer population. If the tracer 
population is a true isolated system of a finite spherical extent 
of r out in equilibrium with the dark halo potential, it follows that 
o~ 2 (r ut) = v(r ou t) = and the boundary term naturally vanishes. A 
similar argument extends to the system of an infinite-extent tracer 
population with a finite-total-mass dark halo, for which r 0M = oo 
also leads to the dropped outer boundary term (then M out = M tot is 
now the 'true' total halo mass). However, if the tracer population 
is pressure-confined and/or the distribution of the observed tracers 
are truncated at a finite outer radius r olll , the outer boundary term 
must remain. In this case, if enough datapoints are available, the 
boundary term may be directly calculated from the observed tracer 
distribution. 



3 EXAMPLES 

We now develop formulae specific to some simple and widely-used 
halo models. 



3.1 The Scale-Free Potential 

The simplest mass model that we consider is (0 < r < r olll ) 



{-V- 



(14) 



Here, a < 1 is the power-index for the scale-free potential, or equiv- 
alently the rotation curve is given by 

2/ \ out 2/ \ 

The central point-mass case is included with a = 1 . If we require 
that the halo density does not increase outwards, then the index is 
restricted to be a > -2, with a = —2 corresponding to a homoge- 
neous sphere. 

For these cases, we have f(r) = 3 - 2/3 + a. If p is further 
assumed to be constant, then / is also constant, and therefore equa- 
tion i ll It reduces to 



(15) 



= (a + 3 - 2/3)— 3g . 

'"out '" oll t 

These are similar to the estimators used bv lWatkins etalJj201Ch 
but the multiplicative coefficient a + y — 2/3 in their equatio n (15) 
is replaced by y => 3. This is because [Watkins et alj j2010j) fur- 
ther assumed a power-law behaviour for the tracer density profile, 
i.e., v oc r~ r , which allowed them to solve for N iot and cr 2 (r). In- 
stead of dropping the boundary term by sending r out — > <x>, which 
is an improper thing to do in a scale-free system, they equated 
Anr'v/Nir) « 3 - y and 07 (7 - 2/3 + a)' l GM/r based on the 
power-law solutions. If we substitute these for the boundary term 
in equation (15) , we recover their equation (15). 

Among the scale-free potential cases, of particular interest are 
the Kepler (a = 1) and the logarithmic (a = 0) potential. For the 
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Kepler potential generated by a central point mass, it is usual to set 
r out = oo and q 2 = 0. Given that /3 is also constant, this leads to 

GM. 



<v 2 r> 



2(2-/3)' 



(16) 



where M. is now the mass of the central point. On the other hand, 
for the logarithmic potential generated by a (truncated) singular 
isothermal sphere, equation U5\ reduces to 

2 i 2\ t 2 

% = <v } - 3g , 

which is actually the same as equation ® with (GM/r) = v 2 being 
constant. 



3.2 The Double-Power Law Halo 

A widely-used family to fit the simulated dark halo density profiles 
is in the form of p(r) <x r a (r% + r i')-^lp . If b > 3, the total 
dark halo mass M tot is finite and jl{r) reduces to the regularized beta 
function with r out = oo. However, for many specific cases the results 
are much simpler. For example, if [p(r)] _1 oc r r (/^+r p ) (3 " r '' p with 
^ y < 3 and p > 0, then 



M(r) 
M,„, 



I _p\-(3-y)//> 



f = 4-2/3- 



3-7 
(r/r )P + 1 ' 



(17) 



That is, the functions p(r) and f(f) are in easily tractable ana- 
lytic form provided that /3(r) is as such. This part icular example 



includes the well-kn own families of the y-s phere l Dehnenl 19931 ; 
Tremaine et alj|l994l) with p = 1 and that of lVeltmannl jl979h and 



Evans & Anl d2005h with y = 2 - p. 

In practice, the simple analytic form of p(r) and f(r) indi- 
cates that the sample mean of pr [ frv 2 is straightforward to cal- 
culate for a fixed set of parameters. Moreover, provided that the 
tracer population extends sufficiently far out (that is, r out s> ro) and 
so Mgat = M(r out ) as M tot , we also argue that it is in general safe to 
drop the boundary term ( formally r oul — » o o, and p. = M/M tot ). 

For instance, for fhe lHernquistl (l990) halo profile, that is, p = 
7 = 1 in equation d 1 Vb . we find that 



GM„ 



(4-2/3 

N 

N 2j 
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1 + 



1 + r/r 
2r, 
r 



(18) 



ro n 



where /3, = /3(r,). If r, <sc r , the observable combination contributes 
like ~ 2(1 - /3i)r 2 l v 2 rj rJ l whereas for r, » ro the contribution is like 
~ 2(2 -yS,)v^r,, which is consistent with the local behaviour of the 
potential at the location of the tracer particle, we also note that the 
mass estimate is dependent upon the choice of the scale length r , 
which is expected. 



3.3 The NFW Halo 

The formal NFW profile dNavarro et alj 19951) , i.e., p~' oc r(r + rf, 
on the other hand possesses infinite total mass, and the proper ap- 
plication of our scheme calls for the truncation of the profile either 
at the radius of the outermost tracer point r out , or at the virial radius 
r v . Let us suppose that the tracers are well-populated so that the 
mass up to the virial radius M v = M(r v ), i.e., the virial mass can be 
effectively estimated. Since M(r) oc m(r/ro) where 

m{x) = ln(l + x) ■ 



for the NFW profile, if we let M out = M v and r out = r v , then 
m(r) 1 / r \ 2 

AW = -H; f(r) = 4-2/3 — — ( — — ) (19) 

where c = r v /ro is the concentration parameter and f = r/ro = 
cr/t\. For a fixed ro, equation i ll It is scaled to 

3c , 



— = (v r h(r)) — g 



(20) 



where h(f) = rf/m(r) whose r-dependence is only via the scaled 
radius f. That is to say, the c-dependence of the result is essentially 
through the overall scale. Here, the boundary term is expected to 
be smalQ and so it is all right to ignore it at the given level of the 
precision. 



4 MASS ESTIMATORS FOR PROJECTED DATA 
4.1 Line-of-Sight Velocity Data 

In many situations, the radial velocities (v r ) - with respect to the 
centre of the halo - of the tracers are not direct observables, but the 
line-of-sight velocities (vf) are. Fortunately, the adjustment of the 
estimator relating to the alternative velocity projections is straight- 
forward. We use the relationship between the line-of-sight velocity 
dispersion (07) and the radial one (07); 

Of = (1 - /3 sin 2 i/>) 07 

where if is the angle between the line of sight towards the tracer and 
the radial position vector of the same tracer from the centre of the 
halo. If the halo is sufficiently far away from us, each line of sight 
towards the individual tracer star or satellite galaxy runs approxi- 
mately parallel to the line of sight towards the halo centre. Then, 
the angle ip is equivalent to the spherical polar angular coordinate 9 
centred on the halo centre. In a spherical system, cr 2 is only depen- 
dent on r and thus by averaging cr 2 over the angle 9 a ip at a fixed 
r, we find that 

-7t/2 



X 71 ' [21 
d^j sin <^ ( 1 - (3 sin 2 i/>) 07 = |l - -yS(r)Jcr 2 . 



Hence the weighted averages of v 2 and v 2 are related to each other 
such that 



3v 2 /i(r) 



(v 2 r h(r)) 



(21) 



3-2/3, 

for any radial weighting function h(r). This is valid for any f3(r) pro- 
vided that the spherical symmetry assumption holds. For the con- 
stant f3 cases, the factor (1 - |/3) is a simple multiplicative constant 
that can be applied after the averaging. 

Finally, the proper form of the mass estimator involving the 
line-of-sight velocities is obtained after substituting equation d21b 
into equation Jilt . 



GM 0U , 13-2/3 + avjr 



3-2/3 p 



'"outs" 



(22) 



1 + x 



where a(r) = 1 — (dln/i/dlnr), which would be a constant a = a 
if p. is given by equation d 14b . 



2 For a virialized system, it is expected that <r 2 (r v ) « 0. In addition, if r v = 
r2oo, then 47ir^v(r v )/W t ot < 4nrlp(r v )/M oM < 3/200 for the tracers that are 
more centrally concentrated than the dark halo. With typical concentration 
parameter values, then 3f 2 c/m(c) < 0.1o" 2 (r v ). 
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4.2 Projected Separation Data 

In more typical cases for an external dark halo, we may only have 
the line-of-sight velocities and the projected distance (R = r sin 6) 
to the halo centre known to within a reasonable precision. Here, we 
would like to find the proper weighting function w{R) of R such that 
(v?w(i?)) = (fr l frv 2 r ), which would replace equation i ll 11 1 with 

GM out = (v 2 w(R)) - 3r ou ,r, 

and lead to the projected mass estimator 



M, 



In practice, the average here would be over the cylindrical region 
with R < R out whereas the one in equation Jilt is over the spherical 
region of r < r m . The distinction is moot if the average is in fact 
over the whole space, i.e., r out = R oM = oo. With a finite cut-off ra- 
dius ^ olll in the tracer population, we proceed by assuming the true 
three-dimensional distribution of tracers is spherically symmetric 
and also cuts off at r out = R out . That is to say, the observed sam- 
ple mean of {v 2 f w(R)) is considered to contain no contribution from 
tracers with r > r ou , and therefore to be a practical estimator for the 
average in the sphere of radius r 0M . 

Provided that the system is spherical and both averages are 
over the sphere, the condition that {vjw(R)) = (v 2 .h(r)) results in a 
integral equation for w(R) at a fixed r; 

r n/2 

rh(r) = r 

(23) 



Jo 



d6 sin 6 (1 - /3 sin 0) w{r sin 9) 
w(R)RdR 



R 2 

1 -P(r)- 



(r 2 -R 2 yi 2 ' 



For h = rf/fl, this is in principle invertible for w(R) if 
lim r _» (dln/i/dlnr) < 3 and (3(r) is finite. We refer the reader to 
Appendix [A] for details. 

In particular, if jS is a finite constant and ft is given by the 
scale-free form in equation d 14b with -2 < a < 1, we find that 

w(R) = T -^- r (24) 



and 



'afl'out 



7t 1/2 r(f + l)a + 3-(a + 2)fi 
4r(2i£) a + 3 -2/3 



where T(x) is the gamma function (the generalized factorial). That 
is to say, the corres ponding mass estimator is in the form of (c.f., 
IWatkins et al]|2O10L eqs. 26 & 27) 



GM„, 



(v)R a ) 



-3 ? 2 



(25) 



Bahc all & Tremaind dl98ll) considered the projected mass estima- 
tor for the central point-mass case. Their results are consistent with 
equation ( 125 b with a = 1 once we drop the boundary term by set- 
ting r ollt = oo, that is, 

32 2-f3 



GM. 



71 4-3/3 
Equation l !25t for a 



y oul J 



(v 2 R). 



on the other hand results in 



(26) 



(27) 



That is to say, if the rotation curve of the spherical halo is flat, the 
line-of-sight velocity dispersion is related to the circular velocity 



(and thus the mass) of the spherical halo such that 3(v 2 ) ~ v 2 ( c.f., 
iLvnden-Bell & Frenkll98iLlEvans et al J 19971: 1 Wolf et alj2010h in- 



dependent of the behavioutj of /?, to an extent that one can ignore 
the boundary term. 

For more complicated mass profiles, an analytic result is in 
general difficult to obtain, except for some special cases. How- 
ever, the special cases do include some interesting examples, one 
of which is the Hernquist halo profile traced by the populations 
with constant f3, for which 



(32 2-/3 R 8 1 -j8r„\ 

h 6 H V 

\ 7t 4 - 3/3 r a 7x2-/3 Rl • 



(28) 



where r is the scale length of the Hernquist halo. The proof is given 
in Appendix|B] together with further examples of mass profiles that 
result in a rational projected mass estimator. 

Strictly speaking, the Hernquist halo cannot be traced by a 
population with j3 > \ at the centre in equilibrium since such a 
population must be cusped at th e centre steeper th an the dark halo 
cusp, which behaves as ~ r _1 dAn & Evansll2006h . However, the 
constant anisotropy needs not extend to the centre and the mass 
estimate with varying anisotropy in the interval (-oo, 1] may be 
understood to be the range of the halo mass consistent with the 
observed line-of-sight velocity data set. 

This projected mass estimator for the Hernquist profile is 
again notably consistent with the scale-free case of equation d25b 
at either extreme, R » ro (a finite total mass; a = 1) or R <k yq 
(r~' cusp; a = -1). This indicates that if the tracers are restricted 
locally in the region where the halo potential can be approximated 
as a power law, equation i25l is a reasonable proxy for the mass 
estimator given proper boundary terms. 

For general mass profiles, the weighting function can be de- 
rived numerically provided that lim r ^o(dln/i/dlnr) < 3 (i.e., a 
cusped halo density profile). While the solution for the general case 
involves a double integral at the least, the function w(R) for a few 
particular cases of constant (3 with an analytic mass profile can be 
obtained through a simple quadrature (see Appendix [A}- For ex- 
ample, the normalized weighting functions W(R) for the NFW pro- 
file - which replaces h{f) in eq. [20] together with the change to 
the line-of-sight veloci ty and R = R/r - for some constant f3 are 
provided in figure 1 of lEvans et ai] feOlCh . For the particular case 
of the NFW profile, it is also possible to derive the power-series 
expansion of W(R) at R = analytically, and also its asymptotic 
behaviour towards R — » oo. In particular, we find at R = that 
W(R) =s ±6 CpR' 1 + $ + 0(R)v/heTe Cp = (1 -/3)/(2-/3) and towards 
R -> oo, that W(R) ~ £/[ln(l +8)-l]~ R/(lnR). 

Finally, the circular orbit model (j3 = -oo) has some special 
points of interest, which are discussed in Appendix lei 



4.2.1 Cored Halo Profiles 

Equation l |25l l is invalid for a = —2 (i.e., an homogeneous sphere of 
radius r out ) because, provided that v 2 4- for R - (i.e., f3 J= -oo), 
the average (v 2 ( IR 2 ) that extends to R = diverges. However, since 
I a= -2ji>-a, also diverges, it actually leads to an indeterminate form. 
In fact, the formal solution for equation < !23t with fi = (r/r out ) 3 
exists in that w(R) = 2(1 - 2/3)r 3 mt 6(R 2 ) where S(x) is the Dirac 
delta - that is to say, letting R" /T(Z + 1) -» i5(fi 2 ) as a -» 2. 



3 Eq. {27) is in fact valid even if/? varies radially. This is because (v 2 ) = 
3<v 2 )for any spherical system independent of the behaviour off}, and there- 
fore eq. (6) is equivalent to eq. (27) for the halo with a flat rotation curve. 
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In practice, the average (v 2 S(R 2 )) is directly related to the tracer- 
number-weighted line-of-sight velocity dispersion along the central 
line of sight, o"l,o, or, 



(v 2 6(R 2 )) = f drvaj = T^X of. 

"tot J "tot 



(29) 



where S is the tracer column density along the central line of sight. 
Therefore, for a = -2 (and /3 ^ -oo), equation \25\ is replaced by 



GM ou 

r oxti 



2(1 -2/?)- 



-3^ 



(30) 



where £ out = MotATf 2 ^,) ' s tne mean column density of the total 
tracer population. 

It can also be shown that the properly derived weight- 
ing function w(R) for any cored halo mass profile with 
lim r ^ (dln/i/dlnr) = 3 contains the Dirac delta. Specifically, if 
lim r ^ i 3 IP- = L is a finite nonzero, then replacing h with h = h - 
2(l-2/3)L/r 2 allows equation i23l to be inverted. If w(R) is the solu- 
tion for this inversion, the final weighting function that reproduces 
(v 2 h(r)) = (v 2 t w(R)) is found to be w(R) = 2(1 - 2/3)LS(R 2 ) + w(R)- 
That is to say, we find that 



GM out 2(1-2j8)L2o< 



(31) 



which differs from the point-mass case in equation J 1 6b by an exact 
factor of two. If the line-of-sight velocity dispersion is used instead, 
we find 

12/2-/? , \ 

in place of equation d22t . 

The calculation for the weighting function suitable for the 
projected separation as a observable is essentially identical to that 
found in Sect. 14.21 as we would like to find the weighting function 
w(R) satisfying {v^w(R)) = ((4 - 2f3)v 2 r), which is to be substi- 
tuted in equation i32\ . This results in the identical integral equation 
d23b with h = 2(2 - yS)r or equivalently fi = 1 . However, because 
of an additional factor of two in equation d32b , the projected mass 
estimator for the self-consistent system with constant [3 is different 
from equation ( 1261 ) again exactly by a factor of two, i.e., 



64 2-yS (v 2 R) 



(34) 



7T 4 - 3/J G 

which encompasses the result of iHeisler. Tremaine & Bahcalil 
dl985h . 



given that 



1 



tL 



■2(1 -2fl)L 



w(R)RdR 
(r 2 -R 2 y 2 



and L = lim r ^ ^Ifi is a finite constant (which further implies that 
lim,_,o f 3 //^ = 2(1 - 2(j)L and thus the above integral equation is 
invertible). 



5 THE SELF-CONSISTENT CASE 

If the tracer density v(r) follows the same functional form as the 
dark halo density p(r) (i.e., the mass-to-light is constant), then fi(r) 
is specified by the integral of v(r) over the volume, and so the prob- 
lem is completely determined. However, for this case the problem 
can be approached through a different simpler route, directly uti- 
lizing the fact that the potential and the tracer density are related 
through the Poisson equation. We find that the self-consistent case 
result in a formally identical mass estimator as the point-mass case 
except for an exact factor of two difference in the associated con- 
stant. 

If v/p is constant, using dM/dr = 4nr 2 p, equation 10 reduces 

to 

dM 2 
~ 

Hence, integrating this on r over [0, oo), we find that 
GM 2 f°° r 4 d , r°° , d (r 4 \ 

where M lot = M(oo) is the total mass. That is to say, 



87tr 4 d , 

(gpo- 2 ). 

G Qdr ' 



dr 



M„ 



-((4-2/?)v 2 r} 



(32) 



The result is valid for an arbitrary functional form for the anisotropy 
parameter /? = fi(r), but if /? is constant, this results in 

GM, nt 



<v 2 r> 



6 DISCUSSION AND CONCLUSIONS 

Here, we have developed the theory of mass estimators. We are mo- 
tivated by instances in astrophysics in which we wish to estimate 
the mass of a dark halo from positions and velocities of tracers such 
as stars, globular clusters, and satellite galaxies. The data sets may 
then be true distances and radial velocities (as for estimating the 
mass of the Milky Way from its satellite galaxies) or may be pro- 
jected distances and line-of-sight velocities (as for the local dwarf 
spheroidal galaxies). In either case, we wish to estimate the mass 
of the dark halo from the kinematics of the tracer population. 

For a given halo density profile, there exists the optimal 
weighting of these kinematic data. We have shown how to find it 
for any given specific density law and different kinds of positional 
and velocity data. This means that the mass within any radius can 
be calculated as the weighted sum of positions and velocities. We 
have worked out the formulae explicitly for a number of important 
cases, including scale-free, Hernquist and NFW haloes. 

Although we have concentrated on general theoretical devel- 
opments in this paper, the performance of the particular scale-free 
estima tors has been tested and was already reported in an earlier 
paper dWatkins et al We have also verified that our estima- 

tors (including the NFW ones) work well against simulation data - 
in which a variety of effects such as halo asphericity, late infall of 
accreted material and lack of virial equilibrium are present. Even 
though these are not taken into account explicitly in o ur estimators, 
none t he less they fare w ell against simulation data jDeason et al.l 
l20ld : lEvans et alJboid) . Finally, applications of our theory to 
estimate the masses of the Milky Way and An dromeda galaxies 
( Watki ns et alj201oT) and the dwarf spheroidals ( lEvans et alj2010h 
are presented elsewhere. 

One notable feature of the mass estimator theories that have 
been developed here is that they do not explicitly depend on any 
properties of the density of the tracer population v. In fact, the mass 
estimator incorporates the information through the definition of the 
average, that is, 



4(2-/3)' 



(u(r,v)) 



— fdV d 3 v f(r, v) u(r, v) = — fd 3 rv(r)u(r) 

JVtot J Mo, J 
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where 



w(r) = —— f d\/(r,v)«(r,v) 
v(r) J 

v(r)= fd^/(r,v); Mot = Jd 3 rd 3 v f(r,v) = fdVv(r) 

and /(r, v) is the phase space distribution function. However, given 
that the sampling of the tracers is statistically random, this formal 
average is estimated by the sample mean, 



<«(r, v)> 



1 N 



That is to say, the effect of the tracer density in our mass estima- 
tor theories is naturally accounted through the spatial frequency 
of sampled tracers, and leaves no explicit dependence on v in the 
consequent formulae. In practice, the choice of the sampled tracers 
may not necessarily be random. If there is some compelling rea- 
son to suspect sampling bias and/or the specific selection function 
is known, the sample mean can be estimated using any additional 
weighting accounting or correcting for the selection bias. However, 
even if a mild selection bias were to be present, the mass estima- 
tor may reasonably be robust without any explicit correction. This 
would be the case if the spatial variation of the quantity to be av- 
eraged for the mass estimator is not strongly correlated with the 
variation of the tracer density itself. 

The conventional way in which mass estimation is performed 
is via the Jeans equations. The procedure is normally as follows: 
first, an assumption is made as to the luminosity density of the 
tracer population; secondly, the data set of discrete velocities is 
binned and smoothed to give the variation of the line-of-sight veloc- 
ity dispersion with radius; thirdly, an assumption as to anisotropy 
is made (often that the anisotropy parameter yS ~ 0) so that the 
line-of-sight velocity dispersion can be converted to the radial ve- 
locity dispersion, and fourthly the spherical Jeans equation is used 
to relate the underlying potential, and hence the enclosed mass, to 
the behaviour of the stellar kinematics. It is worth emphasizing that 
the results obtained for the matter distribution are often not robust, 
as they depend not just on the luminosity profile and the second 
velocity moments, but also their gradients. 

In some sense, the techniques in this paper discard the wealth 
of information contained in the observed data set by taking spatial 
integrals over the whole system. However, there can be a number 
of advantages of this approach, especially if the number of data- 
points are limited. First, it guards against overinterpreting the data 
set which can often happen with use of the Jeans equations. Second, 
the binning of the data, and their subsequent smoothing, are not 
needed. This is actually a great help, as Jeans modelling requires 
derivatives of functions derived from the binned and smoothed data. 
Third, the mass estimators are simple, requiring only weighted 
sums of positions and velocities, as opposed to solution of (at best) 
an ordinary differential equation. 

For these reasons, we expect using mass estimators for dis- 
crete data to be a viable alternative to Jeans modelling. In the limits 
of large numbers of datapoints, we expect the mass estimators and 
Jeans modelling t o yield similarly an swers. This is borne out by 
the calculations of lEvans et al.l ( 1201 Of) for dwarf spheroidals where 
data sets of thousands of radial velocities are available. When the 
the number of tracer datapoints is small, as often happens for esti- 
mating halo masses from satellite galaxies, mass estimators are the 
technique of choice. The precious data need not be smoothed, and 
the estimates of the enclosed mass are robust. We hope our tech- 



niques to be widely used in these instances, as Jeans modelling is 
either too cumbersome or too elaborate. 



ACKNOWLEDGMENTS 

This work originated from the first author's two separate two-week 
visits (February and October 2010) to the IoA (Cambridge), which 
were in part supported by the National Natural Science Founda- 
tion of China (NSFC) Research Fund for International Young Sci- 
entists, as well as the IoA's visitor grant. The authors also thank 
Laura L. Watkins and Alis J. Deason for Monte Carlo tests con- 
cerning the performance of the mass estimators with the Hernquist 
and the NFW haloes. JHA is supported by the Chinese Academy 
of Sciences (CAS) Fellowships for Young International Scientists, 
Grant No.:2009Y2AJ7. 



REFERENCES 

An J. H, Evans N. W., 2006, ApJ, 642, 752 

An J. H, Evans N. W., 2009, ApJ, 701, 1500 

Bahcall J. N„ Tremaine S., 1981, ApJ, 244, 805 

Deason A. J., McCarthy I„ Font A., Evans N. W., Frenk C.S., Be- 

lokurov V., Libeskind N, Theuns T, Gain R., 2010, MNRAS, 

submitted 
Dehnen W„ 1993, MNRAS, 265, 250 
Evans N. W„ An J. H, 2005, MNRAS, 492, 498 
Evans N. W., An J., Deason A. J., 2010, ApJ, submitted 
Evans N. W., Hafner R. M., de Zeeuw P. T, 1997, MNRAS, 286, 

315 

Evans N. W., Wilkinson M. I., Perrett K. M., Bridges T. J., 2003, 
ApJ, 583, 752 

Heisler J., Tremaine S., Bahcall J. N., 1985, ApJ, 298, 8 
Hernquist L., 1990, ApJ, 356, 359 
Kulessa A. S., Lynden-Bell D., 1992, MNRAS, 255, 105 
Limber D. N, Mathews W. G, 1960, ApJ, 132, 286 
Little B., Tremaine S., 1987, ApJ, 320, 493 
Lynden-Bell D., Frenk C. S., 1981, The Observatory, 101, 200 
Navarro J. F, Frenk C. S., White S. D. M., 1995, MNRAS, 275, 
720 

Strigari L. E., Bullock J. S., Kaplinghat M., Simon J. D., Geha M., 

Willman B., Walker M. G, 2008, Nat, 454, 1096 
Tremaine S., Richstone D. O, Byun Y.-L, Dressier A., Faber S. 

M., Grillmair G, Kormendy J., Lauer T. R., 1994, AJ, 107, 634 
Tully R. B., Rizzi L., Dolphin A. E., Karachentsev I. D., 

Karachentseva V. E., Makarov D. I., Makarova L., Sakai S., 

Shaya E. J., 2006, AJ, 132, 729 
Veltmann 0.-1. K., 1979, AZh, 56, 976 (English translation in 

SvA,23, 551) 

Walker M. G, Mateo M., Olszewski E. W., Penarrubia J., Evans 

N. W, Gilmore G, 2009, ApJ, 704, 1274 
Watkins L. L„ Evans N. W., An J. H, 2010, MNRAS, 406, 264 
White S. D. M., 1981, MNRAS, 195, 1037 
Wilkinson M. I., Evans N. W„ 1999, MNRAS, 310, 645 
Wolf J., Martinez G. D., Bullock J. S., Kaplinghat M., Geha M., 

Munoz R. R., Simon J. D., Avedo F. F, 2010, MNRAS, 406, 

1220 



© 2010 RAS, MNRAS 000,[TJ{9] 



8 An & Evans 



APPENDIX A: HOW DO WE INVERT FOR THE 
PROJECTION WEIGHTING FUNCTION? 



r'h ■ 



First, let us rearrange equation d23 b into an equivalent form; 

p 

where 



P Jo 



2 ) 1/2 w(R)RdR + (l -p)G{r) 



G{r) = 



However, we find that 



d7 



r Jo 



{/ -R 2 y n w(R)RdR 



G(r) 



(Al) 



(A2) 



(A3) 



Therefore, differentiating equation iJaTJ after dividing it by r leads 
to a differential equation for G(r); 

/r 3, 



r d7(— )■ 



(A4) 



Given /i(r) and P(r), it is straightforward to solve equation dA4t 
numerically with the boundary condition G(0) = 0. Furthermore, 
equation dA4t can be brought to an exact form 

r<2 d^/r^j 



dr 1 



!dr 



(A5) 



by means of the integrating factor Q(r) satisfying 
dln<2 r d(P\ 



dr 



1 



f dr 



Hence, if Q(r) can be found, it is even possible to bring G(r) to a 
quadrature. Finally, once G(r) is found by some means, the weight- 
ing function w(R) can be obtained though the inverse Abel trans- 
formation of equation l lA2t . i.e., 



w(R) 



2 

Ttfl 2 



r 

Jo 



dr 



dG 



(R 2 - r 2 ) 1 ' 2 dr 

As an example, if /? is a constant, then Q = r s where s 
PI (I -/3)and thus 

G^J-lfdr^ 

1 ~P r Jo dr\ fi I 

The weighting function w(R) is then found to be 

f rdr d 2 r 1 F f(r)r 3+s df 



(A6) 



(A7) 



r R rdr d 2 T 1 C r /(r)r 3+! drl 
J (R 2 - r 2 )!/2 dT 2 J /i(f) J' 



w(7?) - 

V ; 7t(l -/3)tf 2 J (tf 2 - r 2 ) 1 / 2 dr 2 [,- Jo 
For the isotropic case (/3 = 0), this simplifies to 
2 f A rdr ' 

Ttfl 2 



w(R) ■■ 



f 

Jo 



( ^-^ d7[ (3+&) | 



where a = 1 - (dln/i/dlnr). Similarly for p 
weighting function is found to be 



i (i.e., s 



w(R) ■■ 



nR* Jo 



rdr d [ r 

W^7yndrl (2 + a) T 



(A8) 

(A9) 
1), the 

(A10) 



For a system with purely radial orbits (p = 1), equation ( IA5b and 
those derived from it are not valid. However, it is still possible to 
solve for w(R) from equation < |23 1 > or jAlb . which result in 

d 2 



w(R) 



_2_ r R 

nR 2 Jo 



rdr . 

(R 2 - r 2 yn^[ { +a >-~ 



(All) 



The result for the purely circular orbit cases (p = -co) on the other 
hand is obtained by inverting equation iC3\ such that 



w(R) 



A 

TlR 2 



dr 



d if 



(R 2 -r 2 y' 2 dr\fi 



& 



(A 12) 



APPENDIX B: THE PROJECTED MASS ESTIMATOR 
FOR THE HERNQUIST HALO 

Let us think of 



7 = (|— £-_ + 3 + 

■3/3r n2-0 

2 2| 16 2-p r 41-/?r \ 
drr vcr;.\ : — -o 2 + 36 1 H 6 



/ 16 2-p R 4 l-/3r \ 2 

h 3 H vi 

\ tt 4 - 3/3 r n2-0 Rl [ 

An f" 

Jo "" '" r \ 7t 4-3/3 r 



7t2 -p r 



where 



r-n/2 

S„ = dd sin"6»(l -y3sin 2 6»). 
Jo 



We find that 

S 2 = ^(4-3/3); S, = l-\p; S =~(2-p). 

lb 3 4 

Next, 

n 4-3/3 r " n2-p r 



Wd7^ + ^ 



2r r- 2 P dr 



Therefore 



roMot Jo dr J 

• rdrr^iro + rf-lr 2 ^ 2 ). 
Jo dr \ I 



2n 



With the spherical Jeans equation for a constant /3 (eq.[2]wifh Q = 
r 213 ) and the mass profile for the Hernquist halo (eq. [TTJwith y = 
p = 1), we find that 



./ 



2tt 
r N K 



r 



drr 2 vGM tc 



GM, C 
2r 



Similar calculations can also demonstrate the existence of a 
rational projected mass estimator for the particular mass models in 
equation J 17b such that 



GM tot (32 2-p R \ , 
— = \[n4^3p7 +3 h 



for the Jaffe model 

M tot r 



p(r): 



and 
GM„ 

ro 
for 



p(r): 



M(r) 



47T r 2 (r + r) 2 ' M lot r + r 



32 2-p R 8 1 -/3r n 
"n 4 -3/3^ + 7T2-/3 



r o\ 2 



M, ot rg 
2n r(r 2 + r 2 ) 2 



M(r) 
Mm 



ri + r 1 



An analytic example for the cored models that require the Dirac 
delta in the weighting function is found for 



p(r) : 



3M tot r 



M(r) 



An (ro + r) 4 ' Mot W + rj 
whose formal form of the projected mass estimator is given by 
GM totJl 32 2-p R +Q+ 241^^ + 2(1 _ 2 ^ 2) j v? \ 



'0 



7t4-3/3r n2-pR 



© 2010 RAS, MNRAS 000,[TJ{9] 



Tracer Mass Estimators 9 



APPENDIX C: PURELY CIRCULAR ORBITS 

For an extreme scenario, one can imagine that all tracers are in cir- 
cular orbits and their orbital phases and orientations are completely 
random (hence there is no net angular momentum of tracer popu- 
lations in each shell of a fixed radius). This corresponds to the case 
that p = -co everywhere. Since cr 2 = everywhere, the v r -based 
mass estimator is invalid for this case, but the line-of-sight velocity 
based ones are still applicable. Although the final result turns out 
to be the same as the simple limit to/3 — > -oo, we can derive them 
via physically consistent routes. 

Let us start by noting that the line-of-sight velocity dispersion 
of the population in purely circular orbits with random orientations 
is given by 

2 vj . 2 GM(r) . . 
cr. = — sin 6 = sin 6 

e 2 2r 

where v c is the circular speed of the spherical halo at r. First, we 
consider the case that the radial distances r of the individual tracers 
to the halo centre are known. Then, the average of v 2 weighted by 
a function h(r) is found to be 

-71/2 



An 



de r 2 sin 6 vcr 2 h(r) 



Jo 

An r™' , v\h r n/2 
drrv— d6 

N tot Jr m 2 J 



(CI) 



An 

A'toi 



drr v- 



GMh 
~3r~' 



Hence, if one chooses h = 3r/£i where fi{r) is the assumed mass 
profile (eq.U2b, the total mass M out can be isolated by 



M out = - ^ 
G \ y. 



(C2) 



which is consistent with equation J22t in the limit of /? — > -oo. 
The result does not involve the boundary term because a 2 = 
everywhere for the assumed system. 

For the case that only the projected distances R to the halo 
centre are available, we consider the similar weighted average of v 2 
as Sect. 14.21 that is, 

/ 2 \ 4tt r , GM r nl2 , 

(v 2 f w(R))= \drr 2 v de snrW(.R). 

x ' Mot J 2r Jo 

Hence, if w(R) is chosen to satisfy the integral equation 
° 3 dR 



2r 4 _ r r w(R)R 3 
T~X V^R 2 )" 2 



■> G(r), 



(C3) 



the total mass is related to the weighted average via GM out = 
(v 2 w(R)). For the scale-free case (eq. 1 14t. we have w(R) oc R a and 
therefore 



M l): 



4r(^) (v 2 R°) 

'/ 2 r(f + 2) cc 1 ' 



which is valid for a > -A. For a central point mass (a = 1), this 
becomes GM l0l = j^{v 2 R). That is to say, the mass estimate under 
the assumption that the tracers in purely circular orbits is smaller 
by a factor of three and 1.5 respective compared to the case that 
they are in the radial orbits or the isotropic case (see eq.!26t. 

The calculations for the self-consistent case are similar. First, 
we use v/N tot = p/M tot (assuming r iu = and r mt = oo) and 
dM/dr = Anr 2 p to further reduce equation \C\\ to 



1 r° d GhdM 2 
M m Jo 6r dr 



With h = 6r, we find that M to , = 6G~ l (v 2 r), which differs from 
equation iC2l with fi = 1 by a factor of two and is also the limit 
of equation J3 3b as ft — » -oo. With R as an observable instead of 
r, we have G(r) = 4r 4 in place of equation !C3t . which ultimately 
leads to the mass estimator, GM tot = j^(v 2 R), which is the same as 
equation j3A\ for /? — > -oo. 
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